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Abstract 

Analytical formulas are derived to compute the first-order effects pro- 
duced by plane inhomogeneities on the point source seismic response of 
a fluid-filled stratified porous medium. The derivation is achieved by a 
perturbation analysis of the poro-elastic wave equations in the plane- 
wave domain using the Born approximation. This approach yields the 
Prechet derivatives of the P — SV- and SH-wme responses in terms 
of the Green's functions of the unperturbed medium. The accuracy 
and stability of the derived operators are checked by comparing, in the 
time-distance domain, differential seismograms computed from these 
analytical expressions with complete solutions obtained by introducing 
discrete perturbations into the model properties. For vertical and hor- 
izontal point forces, it is found that the Prechet derivative approach 
is remarkably accurate for small and localized perturbations of the 
medium properties which are consistent with the Born approximation 
requirements. Furthermore, the first-order formulation appears to be 
stable at all source-receiver offsets. The porosity, consolidation param- 
eter, solid density and mineral shear modulus emerge as the most sensi- 
tive parameters in forward and inverse modeling problems. Finally, the 
Amplitude- Versus- Angle response of a thin layer shows strong coupling 
effects between several model parameters. 

Kejrwords: Wave propagation, poro-elastic medium, Frechet derivatives 
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I. INTRODUCTION 



The evaluation of the sensitivity of a seismic wave field to small perturbations of the 



material properties is a classical issue of seismology which arises in t 



le solution of forward 



and inverse scattering problems (lAki and Richards 



1980 



Tarantola 



1984 ). Sensitivity op- 



erators are mainly useful for the optimal design of field or laboratory experiments, for the 
interpretation of time-lapse monitoring surveys, and for the development of imaging and 
linearized inversion techniques. In particular, sensitivity operators play a central role in 
least-square inversion schemes using gradient techniques. All these applications, especially 
the latter, call for fast and effective numerical computation methods of the sensitivity op- 
erators. Indeed, the most intuitive approach to compute the perturbational or differential 
seismograms in a structure described by parameters is to use a finite perturbation scheme 
that demands A^ + 1 forward modeling computations. This approach rapidly becomes pro- 
hibitive as the number of parameters increases. The problem can be solved more elegantly 
by first deriving the so-called Frechet derivatives of the seismic wave fields with respect to 
material properties. Although restricted to first-order effects only, this procedure makes 
it possible to efficiently predict small changes of the seismic response resulting from slight 
modifications of the material properties. Furthermore, this solution can be implemented by 
solving only one large forward problem. 



r 



One of the most significant contributions in this area is the work of iTarantolal (Il984j ) who 



applied a first-order perturbation analysis to the elastodynamic wave equation to derive a 
series of general formulas relating the scattered wave field to heterogeneities in an arbitrarily 



complex elastic medium. The same approach was su 
sense fo r layered media in the p l ane-w ave domain by 



jsequently used in a more restrictive 



Pan et al. 



and by 



( 1l988l ) in the acoustic case 



Dietrich and Kormendil (Il990l ) in the P — SV case. The Frechet derivatives ob- 



tained in these different cases are all expressed as combinations of the incident wave field 
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generated by the seismic source at the scatterer location and the Green's functions between 
the scatterer and receiver locations. The structure of these expressions underlines the fact 
that the scattered waves due to perturbations in the material properties (usually density 
and acoustic or elastic parameters) can be interpreted as a wave field generated by a set 
of secondary body forces coincident with the heterogeneities, and determined by complex 
interactions between the incident waves and the medium perturbations. This structure is 
met with all problems of acoustic, seismic or electromagnetic wave propagation in weakly 
inhomogeneous media, and will also be found in more complex situations such as anisotropic 
or poro-elastic media. 

The poro-elastic model, which is the subject of this paper, involves more parameters 
than the visco-elastic case, but on the other hand, the wave velocities, attenuation and dis- 
persion characteristics are computed from the medium's intrinsi c properties without having 



to resort to empirical relationships. Since th e 



( lAuriault et al. 



1985 



Dutta and Ode 



1979 



pioneering wo r 



i of 



Johnson et al. 



Biot 



(119561 ). many authors 



19941 . for example) have intro- 



duced improvements of the poro -elastodynamic equations, either by averaging or by integrat- 



ing techniques. The 



Biot 



(119561 ) theory and its applications is still a field of active research, 
as demonstrated by the large number of current publications devoted to the subject (see, 
e.g., Trifunac, 2006). The forward problem, i.e., the computation of synthetic seismograms 



i n poro-elastic media has b een so 



( ICarciond . 



1996 



Dai et al. 



1995 



ved in different configurations and with several technique s 



Garambois and Dietrich 



2002 



Haartsen and Pride 



19971) 



However, the inverse problem has only been rarely addres sed, and to our knowledge, it has 



never been without first estimating the wave velocities (IBerryman et al. 



2002 



Spikes et al. 



2002 



Chotirod . 



20061 ). Yet, inversion algorithms can provide useful information on the 



material properties, notably permeability and porosity which are the most important pa- 
rameters to characterize porous media. 

The main objective of this work is to extend the methodology used in the elastic case 



( iDietrich and Kormendi 



1990) to obtain the Frechet derivatives for stratified poro-elastic 



media. We consider here a depth-dependent, fluid-saturated porous medium representing 



reservoir rocks or sedimentary layers. The computation of the point source seismic response 
of the layered structure is carried out b y combining the Generalized Reflection a nd Trans- 



missi on Matrix Method (IKennett 



19831 ) with the discret e wavenumber technique 



198ll ). This combination was already implemented by 



Garambois and Dietrich! (120021 ) for 



( Bouchonl . 



the numerical simulation of the co upled s eismi c and electromagnetic wave propagation in 



porous media by using the work of 



Pridd (119941 ) ■ In the following sections, we first present 



the governing equations and constitutive parameters for porous materials before expressing 
the wave propagation equations for depth-dependent media. Next, we develop the analyt- 
ical computation of the Frechet derivatives in the frequency-ray parameter domain for the 
P — SV- and SH-wave cases. Finally, we check the accuracy of the sensitivity operators 
obtained in the time-distance domain, both in an infinite medium and in a complex seis- 
mic model. We conclude with the sensitivity of the seismic waveforms with respect to the 
different model parameters. 



II. WAVE PROPAGATION IN STRATIFIED POROUS MEDIA 



A. Governing equations 



Assuming a e depen dence, where t is time and u the angular frequency, iPridd (119941 . 



20051 ) rewrote Biot's (119561 ) equations of poroelasticity in the form 



V-T = —Uj'^{pU + PfW) 

T = [ KuV -u + CV ■w]I + G[Vu+{Vuf -2/3{V ■ (ul))] 
-p = CV-u + MV- w 



-VP = —Uj'^PfU — UJ^pW 



(1) 



where u and w are the solid average displacement and the relative fiuid-to- solid displacement, 
respectively. More precisely, defining Ug and Uf as the displacements of the solid and fluid 
phases of a porous continuum, we can write u Ug and w = 0(Us — Uf), where is the 
porosity. P represents the interstitial pressure and r is the 3x3 stress tensor, p denotes the 



density of the porous medium which is related to the fluid density p/, solid density ps and 
porosity via the relationship 



P 



1 



Ps + (ppf ■ 



(2) 



The undrained bulk modulus Ku is defined under the condition w = 0. G is the shear 
modulus of a dra ined or an undrained medium as it is independent of the fluid properties 
). The fluid storage coefficient M represents the amount of fluid a sample 



( IGassmann 



1951 



can accumulate at constant sample volume. Blot's C modulus is a mechanical parameter 
describing the variation of the fluid pressure due to a change of the sample volume in an 
undrained medium. At low frequencies, these parameters as well as the Lame parameter 
Xu defined below are real, frequency-independent and can be expressed in terms of the 
drained modulus -ft^n, porosity 0, mineral modulus of the grains Kg and fluid modulus Kf 



( IGassmann 



195l|): 



Ku 
C 

with A 



1-(1 + 



K 



D 



Kf 



K 



D 



K. 



(1 + A) 



K 



u 



G 



:i+A) 



1-6 Kf 



K. 



1 - 



M 
K 



Kf 



;i+A) 



(3) 



D 



K. 



It is also possible to link the bulk properties K^, and G to the porosity and con- 



(Bemer et ai. 


2004; 


Pride, 


2005) 



Kn = K, 



1 - 



and 



G — Gfi 



1 - 



(4) 



1 + c, " " ^ 1 + 3 0/2 ■ 

Equations (jlj) have the merit of being very simple and introduce only two additional param- 
eters, namely, the shear modulus of the grains Gs and the consolidation parameter Cg. The 
latter mainly depends on the cementing properties of the grains, but also on the pore shape. 
The value of the consolidation parameter is typically between 2 to 20 in a consolidated 
medium, and can be very much greater than 20 in an unconsolidated soil. 



Finally, the wave attenuation is expla ined by Darcy's law which uses a complex, 
frequency-dependent dynamic permeability (jjohnson et all Il994j ) : 



uj k{(jj) 



with k{ijj) = fco/ 



(5) 



The dynamic permeability kiuS) tends toward the hydrogeological (dc) permeability k^ at 
low frequencies where viscous losses are dominant. It includes a correction accounting for 
the inertial effects at higher frequencies. These two domains are separated by the relaxation 
frequency 



(6) 



where t] is the viscosity of the fluid. Archie's law ( F = ) expresses the electrical 
formation factor F in terms of the porosity (p and cementing exponent m whose value is 
between 1 to 2 depending on the pore topology. Parameter rij is conside red co n stant and 



equal to 8 to simplify the equations. We refer the reader to the work of 
more information on the parameters used in this study. 



Pridel fl2005f ) for 



B. Coupled second-order equations for plane waves 



The horizontally layered model lends itself to a number of analytical developments if one 
performs a plane wave decomposition of the wave fields represented by equations ([1]). This 
procedure in volves a s e ries o f changes of variables and integral transforms which are described 
in detail by 



Kennett 



new variables (Hudson 



1983) in the elastic case. As in the elastic case, the introduction of 



19691 ) leads to a useful decomposition into P — SV- and SH-wave 
systems in cylindrical coordinates. By applying the whole sequence of transformations and 
arranging the terms, we find that the governing equations of the P — SV-wave system in 



depth-dependent poro-elastic media reduce to the following system of second-order equations 



in the angular frequency uj and ray parameter j9 domain: 



d_ 

dz 



dU dW 
{Xu + 2G) — + C^-ujp{XuV + CX) 
oz oz 



dV 

-ujpG — 
oz 



+ uj^ [pU -p^GU + pfW 
d 



J- Ir 
F2r 



dz 



dV 

G — + ujpGU 
oz 



+ ujp 



oz Oz 



+ [ pV - p^ {\u + 2G)V - p'^ C X 



d_ 

dz 

ujp 



dU dW 

C—-ujpCV + M- up MX 

oz oz 



+ [pfU + pW 



dz dz 



+ [pfV + pX -p'^{CV + MX) 



(7) 



In these equations, U = U{zji,uj; zs) and V = V{zr,u;zs) respectively denote the verti- 
cal and radial components of the solid displacements. Similarly, W = W{zr,uj; zs) and 
X = X{zfi,u!; zs) respectively denote the vertical and radial components of the relative 
fluid-to-solid displacements. Variables zr and zs stand for the receiver and seismic source 
depths. Equations ([7]) are valid in the presence of body forces Fi = [Fiz{zs, u), Fir{zs, uj)]'^ 
and F2 = [F2z{zs,uj), F2r{zs,uj)]'^ defined by their vertical {z index) and radial (r index) 
components: force Fi is applied on an average volume of porous medium and represents a 
stress discontinuity, while force F2 is derived from the pressure gradient in the fluid. 

We can then cast equations ([7]) in the form of a matrix differential equation as 



PSV PSV 
Li (4 



where Q 



U 
V 

w 

X 



and F 



Fiz 

Fir 
F2, 
F2r 



(8) 



PSV 



is a differential operator given by 



j^psv d 

dz 



, ^PSV d , ^PSV 



-iUp [M2 



PSV ^rp d 2 ' - -PS"^ 



-+.^(M3 



p'mr) (9) 



8 



PSV 

where the M ■ , i 



1..4 are 4x4 matrices defined by 



Ml 



Xu + 2G C 






G 








c 





M 


















p Pf 

p pf 

Pf p 

Pf p 



-Xu -C 
G 

-C -M 


a 

Xu + 2G C 

C M 



(10) 



Apart from the dimensions of the matrices, we may note that equations ([9]) and ffTOl) are 
very similar to the expressions obtained in the elastic case. 

Formally, equation ([8]) admits an integral solution for the displacement fields in terms 
of the Green's functions of the problem. For example, the vertical displacement U at depth 
zr and frequency u is given by 



U{zn,uj)= / [G\l{zji,uj;z')F,,{z',uj)+Gil{zji,uj;z')F,,{z',uj) + 
Jm 

G\l{zn, oo; z') F^r{z' , u) + GjUzn, u; z') F^riz', oo) ] dz' 



(11) 



where G^!j{zji,uj, zs) is the Green's function corresponding to the displacement at depth 
Zfi of phase i (the values i = 1,2 correspond to solid and relative fiuid-to-solid motions, 
respectively) in direction j {z or r) generated by a harmonic point force Fki{zs, uj) {k = 1, 2) 
at depth zs in direction I {z or r). A total of 16 different Green's functions are needed to 
express the 4 displacements U, V, W and X in the P — SV-wave system (4 displacements 
X 4 forces). 

The integrals of equation (fTTj) are taken over the depths z' of a region A4 including the 
body forces Fi and F2. In the case of a vertical point force at depth zs, the expressions of 
the forces become 

Fi,{zs,uj) = 6{z - Zs) Si{uj) \ F2,{zs,uj) = 6{z - zs) S2{uj) 

and < 

Fir{Zs,Uj) = \ F2r{Zs,Uj) = , 



[12) 



where 5*1(0;) and 5*2 (cj) are the Fourier transforms of the source time functions associated 
with fo rces and F2. Assuming tha t the amphtudes of both forces are similar (||Fi|| ~ 



F2II) (iGarambois and Dietrich 



2OO2I ). we take S{iu) = Si{iu) = 5*2(0;). The displacement 
fields for a vertical point force can then be written in simple forms with the Green's functions: 

U{zr,uj;zs) = [G\1{zr,uj;zs) + Gjl{zR,uj; zs) ] S{uj) 
V{zr,u;zs) = [GY^{zr,uj;zs) + Gf^{zR,uj; zs) ] S{uj) 

(13) 

W{zr,uj;zs) = [Gll{zR,uj;zs) + GII{zr,uj; zs) ] S{uj) 
X{zr,uj,Zs) = [GI^^{zr,uj,Zs) + GI^^{zr,u;,Zs)] S{uj) . 
The displacement fields corresponding to a horizontal point force and to an explosive point 
source are similarly defined. Explosions would be represented by Green's functions G^^^ 
{i = 1, 2; j = z,r; k = 1, 2) representing the radiation of an explosive point source E. 

The 5"-^ case is treated in exactly the same way as the P — SV case. The corresponding 
second-order differential equations of motion are 



d_ 

dz 



G 



dT 



dz 



+ uj^[ -p^GT + pT + pfY] 



(14) 



Fu 

c 

F2t = uj^PfT + pY], 
where T = T{zr, u; zs) and Y = Y{zr, u; zs) stand for the transverse solid and fiuid-to-solid 
displacements; Fu{zs,uj) and F2t{zs,Lj) are the transverse components of body forces Fi 
and F2. 

As before, we can rewrite equations f|T^ in matrix form as 



j^SH^SH ^ ^ ^^^^^ ^SH 



T 
Y 



and F" 



Fit 

F2t 



(15) 



SH 



Here, L is a linear operator defined as 



_ SH d I ^ ^SH d \ r, SH 9 , ^SH 

where the M^",2 = 1..4, are 2x2 matrices defined by 



(16) 



G 





-» «-SH 


P 


Pf 






, M2 = 












Pf 


P 



(17) 
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III. FRECHET DERIVATIVES OF THE PLANE WAVE REFLECTIVITY 



A. Statement of the problem 

The Frechet derivatives are usually introduced by considering the forward problem of 
the wave propagation, in which a set of synthetic seis mograms d is co mputed for an earth 



model m using the non-linear relationship d = f(m). iTarantolal (119841 ) uses a Taylor series 
expansion to relate a small perturbation 6m in the model parameters to a small perturbation 
Si in the wave field 

f(m + 5m) = f(m) +D (5m + o( ||5mf ) or 5f = D 5m (18) 

where T) = di / dm is the matrix of Frechet derivatives. 

Our aim is to compute the various Frechet derivatives corresponding to slight modifi- 
cations of the model parameters at a given depth. Considering for instance the density at 
depth z, this problem reduces, in the P — SV case, to finding analytical expressions for the 
quantities 

.PSV, I . dU{zR,uj;zs) 

Ai {zr,uj;zs\z) 

Al''^{zn,uj;zs\z) 

m 

izR,uj;zs\z) 
AI^^{zr,uj; zs\z) 



dp{z) 




dV{zR,uj] 




dp{z) 




dW{zR,uj 


;^5) 


dp{z) 




dX{zR,uj; 





dp{z) 

We can similarly define the Frechet derivatives C^^^ , T^T^ i ^i^^ i ^i^^ ^iid GJ^^ , 

i = 1..4, for model parameters p/, p, C, M, Xu and G. This is the natural choice of 
parameters to carry out a perturbation analysis because of the linear dependence of these 

^ PSV 

parameters with the wave equations. We also introduce the set of Frechet derivatives Aj^ , 
B^^^ , Cf^^, V^^^ , Sj^^ , ^J^^ , GJ^^ and H^^^ , i = 1..4, corresponding to model parameters 
Ps, Pf, ko, 0, Ks, Kf, Gs and which we will use in a second stage, and which are more 
convenient to use as physical parameters of the problem. The sens i tivity operators are 



derived by following the procedure presented in 
case. 

11 



Dietrich and Kormendil (119901 ) for the elastic 



B. Perturbation analysis 

We first present, with some detail, the perturbation analysis for the P — SV case before 
addressing the simpler SH case. We consider small changes in the model parameters at a 
given depth z that result in small perturbations AQ = [5U , 5V , 5W , 5XY of the seismic 

PSV 

wave field and in a modified seismic response Q' = [[/', V' , W, X']'^. The latter can 
be written as 

, PSV pgv PSV 

Q' = Q + AQ (20) 

PSV 

by assuming that the magnitudes of the scattered waves AQ are much smaller than those 

PSV 

of the primary waves Q . Considering for instance the first component of the above vector 
equation 

U'{zr, u; zs) = U {zr, uj] zs) + 5U {zr, u; zs) , (21) 
we can write the scattered displacement 5U as 



5U{zr,uj;zs) = / {zR,uj;zs\z)6p{z) + {zr,uj; zs\z) 6pf{z) + 

Jm 

PSV PSV 

{zR,uj;zs\z)5p{z) {zr,uj;zs\z)5C{z) + 

PSV PSV 

£^ {zr,u;;zs\z) 5M{z) + {zr,uj; zs\z) 5Xu{z) + 
gr''{zR,u;;zs\z)6G{z)]dz . (22) 

With the model parametrization adopted (i.e., a linear dependence of the model parameters 
with the wave equation), the perturbation analysis can mainly be done in symbolic form. 

PSV 

Indeed, the wave operator L in the perturbed medium can be written as 

L'''^ = L^^^ + AL^^\ (23) 



so that equation ([8]) becomes (IHudson and Heritagel . Il98ll ). 



/_ PSV * T PSV\ /„PSV * ^PSV\ -r-,PSV 

f L + AL j (Q + AQ j = F (24) 

We then use the single scattering (or Born) approximation to solve the above equation for 
AQ''^^ under the assumption (already used above) that ||AQ^^^|| ^ ||Q''^^||. Combining 

12 



this approximation with equation ([8]), we obtain 



, PSV „PSV 

L AQ ~ 



A 1 ^PSV 

AL Q 



AF 



PSV 



(25) 



This matrix equation shows that the scattered waves AQ due to perturbations of the ma- 

PSV 

terial properties can be interpreted as a wave field generated by secondary body forces AF 

PSV PSV 

defined by the interaction of the incident waves with the heterogeneities (term AL Q ). 
Moreover, this wave field propagates in the unperturbed medium represented by the wave 
operator L''^^. Consequently, equation (125|) has for each of its components a solution similar 

PSV PSV PSV PSV 

to equation (fTTj) . by substituting AF for F and AQ for Q . Considering again 
the scattered displacement 6U, we have 



6U{zr,uj]Zs) 



G\l{zR, uo- z') 6Fi,{z', Lu) + GjlizR, u; z') 6F2,{z', u;) + 
M (26) 

G\:izR,uj;z')SF,riz',cu) + Gllizn,co;z')5F2riz',u) ] dz' , 



where the secondary Born sources ^-Fi^, 6 Fir, 5-^22, SF2r are obtained from equation ([7]): 



AF 



^SFi\ 

SFir 
\ SF2r ) 



( _d_ 

dz 



dU dW 
{6\u + 26G) — + 6G^-ujp {6\u V + 6G X) 
oz oz 



dV 



+ uop5G—-uo'^ [ 5pU -p^dGU + 5pfW] 



d_ 

dz 

d_ 

dz 



-cop 



dz 
dV 

SG— + ujpSGU 
dz 



cop 



dz 



dW 

dz 



co^ [6pV-p^ {6Xu + 2 6G)V -pHGX 



dU 

SG [ —-copV ] +SM 



co^ [ 6pfU + 6pW ] 

,^dU ,,,dW 
6G— + 6M—- 
dz dz 



dW 

dz 



copX 



CO 



\ 



[ 6pfV + 6pX -p"^ (SGV + SMX) 



J 



(27) 



By inserting these expressions into equation (!26l) and integrating by parts in order to separate 
the contributions in 6p, 6pf, 6p, 6G, 6M, 6Xu and 6G, we obtain an integral representation of 
the scattered wave field 6U that we can directly identify with equation (j22|) to get the Frechet 
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derivatives , , , , , J^^ 

Ar = -u;'[UGll + VGl:] 

-uj^ [WG\l + UGll + V Gf, ] 

-u'[WG\l + XG\:] 

— ujpX 

— ojpV 



,PSV ^PSV 



and Qi corresponding to displacement U : 



^1 



dW 
dz 
dW 



Iz 



dz 
dU 
dz 



dV 

— + ujpU 
oz 



dG\ 

dz 
dGll 

dz 
dG\l 
dz 
dG^ 
dz 



-ujpGi: 
-upGi: 
-upGi: 



+ 



dU 
dz 



— ujpV 



2z 



dGj 

dz 



-upG{ 



2r 



(28) 



+ ujpGi 



Iz 



dU dGll 2 2 T r 



Here, U = U{z,u; zs), V = V{z,u; zs), W = W{z,lj; zs) and X = X{z,u;zs) denote 
the incident wave fields at the model perturbations. These wave fields can be expressed in 
terms of Green's functions using equations (|T3|) for a vertical point force, and similarly for 



a horizontal point force or for an explosive point source. Expressions Gfj = Gfj{zR,u!; z) 
represent the Green's functions conveying the scattered wave fields from the inhomogeneities 
to the receivers, as noted before. A total of 32 Green's functions (16 for up-going waves and 
another 16 for down-going waves) are to be computed to completely solve these equations. 

The Frechet derivatives for the radial displacement V are easily deduced from expressions 
(125]) by changing G^f, to Gf^.. In the same way, the Frechet derivatives corresponding to the 
vertical and radial relative fiuid-to-solid displacements are obtained by substituting G2j for 
Gij. As a first verification of the formulas, we note that the expressions corresponding 
to perturbations of parameters p, Xu and G {A^^^ , J-'J^^ and QJ^^ , respectively) are in 



perfect agreement with the Frechet derivati ves formulas computed 



Lame parameters A and /i in the elastic case (jPietrich and Kormendi 



' or the density p and 



1990). 



C. Frechet derivatives for relevant parameters 



As mentioned above, the set of 7 parameters (p, pj, p, C, M, Xu and G) used in 
the perturbation analysis was primarily chosen to considerably simplify the derivation of 

14 



expressions ( !28|) . In practice, and as with the Lame parameters in the elastic case, it is 
more convenient to consider model parameters which are easier to measure or estimate. We 
introduce here a new set of 8 parameters, namely p^, pj, k^, 0, Kf, Kg, Gg and which 
are more naturally related to the solid and fluid phases. Furthermore, these parameters are 
independent from each other in terms of their mechanical or hydrological meaning. They 
are all real and frequency independent, contrary to p. Fluid viscosity is not considered here 
as the Frechet derivatives with respect to the fluid viscosity and to the permeability are 
strictly proportional. We concentrate on the permeability by stressing that the comments 
made for the permeability will also be valid for the fluid viscosity. 

To obtain the Frechet derivatives corresponding to the new set of parameters from 
expressions (128|1 . we construct the 8x7 Jacobian matrix J whose coefficients are formally 
defined by 

where = 1..7, stands for one of the parameters p, p/, p, C, M, Xu or G, and where 
p'j,j = 1..8, represents one of the parameters ps, p/, ko, cj), Kg, Kf, Gg or Cg. The elements 
of matrix J are determined from equations (Ej) to ([6]). 

The Frechet derivatives with respect to the new set of parameters for the vertical component 
of the solid displacement U are then defined by 



PSV 



;psv 



PSV 



/^PSV 

PSV 



(30) 



\n: j 

The combinations involved in equation ( l30ll significantly complicate the expressions of the 
Frechet derivatives corresponding to the new set of parameters. For sake of clarity and 
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simplification, we introduce the following quantities 



— A rij — Ai — 

V 

a = (1 + c,0) Ks + Cs{l - 0) Kf 

6 = 02 {4 (2 + 3 Cs) - 3 (2 + 3 Q 0)^ K, (1 + c) [i^. 



(31) 



With these parameters, the final expressions of the Prechet derivatives Ji^^, B^^^, C^^^, 



Tfi^^ Sl^^ ^ J^i^^i Qi^^ and with respect to model parameters p^, p/, /cq; 0; -^s; 



Gg and are 



^1 = (1 - 0) ^1 

^PSV , .PSV , „PSV , 2 + VL PSV 



PSV 



i n J + 2 



UJr. 



V 



PSV 



I 02 {p, - pf) (2 + 3 c, 0)2 aT - 2 a2 02 (2 + 3 c,) g"^^ 



a? 02 (2 + 3 c, 0)2 

- (2 + 3c,0)2 [(l + C,0)2X, + C,(l-20-C,02)ir; 



PSV 



- 02c,(l + c,)(2 + 3c,0)2(K,-i^y)i^,i^y + -hT^- 
sT =^-^{^-s{l + c.) K,- VT + c. (1 + c. 0) i^/ 



^0 
+ 



(32) 



(1 + C,0)ii:,2 + C,(c, + 0)X/ + 2C,(1-0)X,X;] TI^^ } 



j^r = 47 1 (i+c.0)ii's'f0(i+c.)pr+ (i+c.0)^r 

(c, + 0) (1 + c, 0) - c, (1 - 0)2 ] X,2 } 



^1 



PSV 



= 2 



1-0 
2 + 3 
1-0 



2 ^PSV ^PJ 



a2 0(2 + 3c^0)2 
— c — oa (p Gs Gi > 



( of PSV P' 

I (2 + 3c,0)2i^,i^^ ^0(i^,-ir^) -i^/^i' 



We note that the Prechet derivatives and with respect to fiuid density p/ 

and permeability /cq a-re complex due to the role of these parameters in the attenuation and 
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dispersion of seismic waves. 

In addition, formulas ([32D can be further simplified if source and receivers are located at the 
same depth zq = zr = zs- In this case, we can take advantage of the representation of the 
incident wave fields f/, V, W and X in terms of the Green's functions (see equations [T3|) . 
and use the reciprocity theorem: 

G\]{zr,uj-zs) = G%{zs,u;-zr) . (33) 

The number of Green's functions required to describe the wave propagation then reduces 
from 32 to 16. These simplifications are straightforward and are not developed here. 



D. SH case 



We follow the same procedure as above to derive the Frechet derivatives of the solid 
displacement T and relative fiuid-to-solid displacement Y in the SH case. We denote by 
jt^ ., B^^ , Ci^, V^^, , J^^^ and Q^^, respectively, the Frechet derivatives with respect 
to model parameters p, pf, p, M, Ku, G and G, where subscript i = 1 refers to the solid 
displacement T and subscript i = 2 refers to the relative fiuid-to-solid displacement Y. We 
also introduce the notations Af^ , C^^, , , , ^™ and T-Cl^ for the Frechet 
derivatives relative to our alternative set of model parameters ps, pj, ko, cj), Ks, Kf, Gg and 
Cs- The perturbation analysis leads to the following expressions of the Frechet derivatives of 
the transverse solid displacement T: 

Br = -u;'[YG\i + TGli] 
Cr = -uj^YGll 

= (34) 

s!"" = 



^'"^ = 



dT dG]' 



SH U\^lt 2 Jlrp ^It 
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We note, as before, that the displacements T and Y can be expressed in terms of the 
Green's functions G^*, where subscripts i = 1 and j = 1 relate to solid displacement and 
force, whereas subscripts i = 2 and j = 2 relate to relative fiuid-to-solid displacement and 
force. Subscript t stands for the tangential displacement or force. The transformation of 
these expressions with the Jacobian matrix J defined in equation (12^ finally yields: 

AT = (!-</') Asm 

-SH 2 + Q 

B- =<pAsHi + Bshi + —^FCshi 



(mj + 2 — )?7 



f>T = (Ps- Pf) Ashi - 2 Gs . \ Gshi 



{2 + ?>Cs(t>y (35) 



S^"" = 



^SH _ 1- (j) 

'^^ = -^^'^{2 + 3cs<py^''' 
We note that the Frechet derivatives with respect to Kg and Kf are zero since these param- 
eters have no influence on shear waves. 



IV. NUMERICAL SIMULATIONS AND ACCURACY TESTS 
A. Frechet derivatives vs discrete perturbations 

As mentioned in the introduction, the Green's functions for layered media are computed 



with the Generalized Reflection and Transmission Matrix Method of 



Kennett and Kerry 



(Il979l ) which yields the plane-wave response in the frequency-ray parameter (or horizontal 
wavenumber) domain. In our numerical applications, the Frechet derivatives are first calcu- 
lated in the frequency-wavenumber domain before being tra nsformed into t he time-distance 
domain with the discrete wavenumber integration method (IBouchoJ . 1981 ). 

In order to test our analytical formulations and assess their limitations, we compare 
the differential seismograms computed with the first-order Frechet derivative approach with 
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TABLE I. Model parameters of the infinite medium used for the numerical tests of the 
Frechet derivatives. 



( ) A;o (m2) pf (kg/m^) (kg/m^) (GPa) (GPa) Kf (GPa) c, ( ) 
0.20 10-12 xoOO 2700 35 25 2.2 50 

the seismograms obtained by introducing discrete perturbations in the medium properties. 

Thus, considering for instance the vertical component of the solid displacement in the P—SV 

case, the partial derivative dU/dpj with respect to parameter pj (j = 1..8, according to the 

parameter set considered) can be approximated by the following finite difference expression: 

dU{zR,u;;zs) ^ [/p.+Ap. (^^^ ^. _ up, (^^^ ^. 
dpj Apj 

where Apj represents a small perturbation of parameter pj. 

The similarity between the seismograms computed with the two approaches indicated in 

the left- and right-hand sides of equation (!36ll is evaluated from the correlation coefficients 

between the traces. 



B. Uniform medium 

We first consider the simple case of small perturbations Sps, Spf, Sk^, Scf), SK^, SKf, SGs 
and 6cs within a thin slab embedded in an infinite uniform medium. The slab thickness is 
of the order of one twentieth of the dominant wavelength of the P- waves (i.e., 1 m). The 
amplitude of the relative perturbations Apj/pj is 10 % for each of the parameters considered. 
Source and receivers are located at the same depth, 50 m above the model perturbation. The 
parameters of the uniform model are listed in table 1. The seismic response of the thin slab 
is considered as the reference for the comparisons with the Frechet derivative seismograms. 

The simulations shown in figure [1] include Pfast-, Psiow- and S'-waves whose computed 
velocities are respectively equal to 2250, 130 and 750 m/s at a frequency of 85 Hz. The slow 
P-waves are not visible, but three refiected waves (compressional PP, converted PS and 
SP, and shear SS) are easily identified in the four sections displayed in figure [TJ It is seen 
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that a small perturbation of the fluid modulus Kj has no influence on shear waves. The 
same behavior is observed for the solid modulus Kg (not shown). On the contrary, slight 
changes in the other parameters mainly generate SS reflections, as noted in particular for 
the porosity (p, mineral shear modulus Gs and permeability ko in flgure [TJ We also observed 
that the differential seismograms are very similar for the following pairs of perturbations: i) 
consolidation parameter and porosity; ii) fluid and solid moduli; in) fluid and solid densities. 
In addition, we found that the correlation coefficients between the Frechet derivative and 
discrete perturbation seismograms are greater than 99 % for most model parameters at all 
source-receiver offsets. The only exception concerns the Frechet derivative with respect to 
the permeability which shows correlation coefficients between 60 and 95 % in the P — SV 
case, depending on the source-receiver offset. However, this operator appears more accurate 
in the SH case. By and large, the tests performed in a uniform medium validate our 
analytical expressions derived in sections III.C and III.D. 

We further check the accuracy and stability of the first-order sensitivity operators by 
modifying the amplitude of the discrete perturbations, with the following results: i) the 
Frechet derivatives with respect to parameters that only influence P-waves are more accurate 
than the Frechet derivatives with respect to parameters that influence both P- and S- 
waves. ii) Strong perturbations of the solid and fluid moduli Kg and Kf do not produce 
any distortion of the waveforms, but merely result in a global increase of the amplitudes of 
the discrete perturbation seismograms. Hi) For strong amplitude perturbations, the Frechet 
derivatives are more stable at near offsets (i.e., for small angles of incidence) than at large 
offsets. We interpret this observation as being due to the nonlinearity inherent to large 
offsets where the wave flelds interact more strongly with the subsurface structure. However, 
as an exception to this rule, the Frechet derivatives with respect to permeability appear more 
stable at large source-receiver offsets. This may be explained by the fact that a perturbation 
of the permeability mainly influences the wave attenuation and dispersion and therefore has 
a stronger effect for longer travel paths. The deterioration of the seismograms correlation 
with decreasing ko (and corresponding decrease in wave attenuation), is another indication of 
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the weak influence of tlie permeability on the seismograms. iv) For strong perturbations, the 
accuracy of the Frechet derivatives deteriorates for specific offsets corresponding to critical 
angles. This reduced accuracy manifests itself by modifications of the relative amplitudes 
of P— and 5*— waves rather than by waveform changes, v) When checked against discrete 
perturbations of positive and negative amplitude of the same magnitude, the first-order 
approximations do not show exactly the same accuracy. In general, the Frechet derivative 
seismograms obtained for positive perturbations display a better accuracy, vi) As a general 
rule, the first-order approximations appear remarkably accurate for amplitude perturbations 
up to 20 % in absolute value. 

We now consider the robustness of the Frechet derivative seismograms with respect to 

Distance (m) Distance (m) 

50 100 150 200 50 100 150 200 



i 




Frechet deriv. Gs modulus Frechet deriv. Permeability 



FIG. 1. Seismic sections corresponding to the first-order perturbation of the vertical dis- 
placement with respect to porosity 0, fluid modulus Kf, mineral shear modulus Gs and 
permeability ko in the uniform and infinite medium described in table 1. The seismic exci- 
tation is a vertical point force. 
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the thickness of the perturbed layer. In our uniform model, the wavelengths Xp and Ag 
corresponding to P- and S-waves are respectively equal to 26 m and 9 m at the domi- 
nant frequency of the Ricker wavelet used in the simulation. Our computations show that 
the Frechet derivative seismograms are very well correlated with the discrete perturbation 
seismograms until the thickness of the perturbed layer reaches about 20% of the dominant 
wavelength Xp (that is, ~ 5 m) for parameters Ks and Kf, and 20% of the dominant wave- 
length Xs (that is, ~ 2 m) for all other parameters that influence both P- and S'-waves. Thus, 
we observe that the Frechet derivatives with respect to Ks and Kf are more robust than the 
other expressions with respect to departures from the "small and localized perturbation" 
assumption of the Born approximation. 

C. Complex model 

We used the 16-layer model shown in figure [2] to numerically check the stability and 
accuracy of the Frechet derivative formulas in a more complex structure. In this model, the 
perturbed layer is at a depth of 50 m, source and receivers being located near the surface. 
Figure [3] presents the seismic sections obtained with the Frechet derivative and discrete 
perturbation methods for slight changes of the solid density ps and permeability ko. On the 
whole, the comparison of the waveforms obtained with both methods is very satisfactory in 
spite of some differences observed at small offsets for the perturbation of the permeability. 

Figure H] shows that the Frechet derivative seismograms are remarkably accurate as 
long as the perturbation amplitude remains small. In this case, the maximum acceptable 
perturbation amplitude is approximatively 10 % of the model parameter value. When the 
perturbation amplitude is increased beyond this limit, the P— waveforms remain practically 
unchanged whereas the 5— waveforms are distorted. In all cases, no variations in travel 
times are observed. 

We also checked the behavior of the Frechet derivative operators relative to the thickness 
of the perturbation layer, or equivalently, relative to the central frequency of the wavelet 
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used in the computations. As for the uniform medium investigated before, we note that the 
first-order approximations remain very accurate as long as the layer thickness layer is lower 
than Ap/5 or A5/5 depending on the model parameter considered. 




0.2 0.4 10 20 1000 2000 10 20 

( ) (10-^^ m^) (kgW) (GPa) 



FIG. 2. 16-layer model used for the numerical tests of the Frechet derivative formulation in 
a stratified medium. The figure shows the distributions of the eight model parameters as a 
function of depth. 

V. SENSITIVITY STUDY 

A. Relative influence of the model parameters 

In this section, we assess the relative influence of small modifications of the model 
parameters on the different components of the seismic wave field. For this, we determine the 
maximum amplitude of the seismograms computed with the discrete perturbation method 
for all pairs of vertical and horizontal forces and displacements. These computations are 
done with the model parameters of table 1, aside from the consolidation parameter Cg which 
is given a value of 20 in figure [5] and 100 in figure [61 Perturbation depth, source and 
receiver locations and maximum offset are identical to those used in section IV. B. The 
amplitudes thus obtained are multiplied by the parameter variation (Ap) to obtain the 
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displacement change Af/ of equation (l36i) . and are normalized with respect to the maximum 
value found. The same quantities were computed with the Frechet derivative approach to 



500 




(a) Frechet deriv. solid dens. 



Distance (m 
200 300 



(b) Discrete pert, solid dens. 



500 




(c) Frechet deriv. perm. 



(d) Discrete pert. perm. 



FIG. 3. (a) and (h) Seismic sections obtained with the Frechet derivatives and with the 
discrete perturbation methods for a perturbation of the solid density at z = 50 m in the 
16-layer model depicted in figure [2l The seismograms represent the vertical displacement 
generated by a vertical force (P — SV wave system), (c) and (d) Same for a perturbation 
of the permeability k^. In this case, the seismograms represent the horizontal transverse 
displacement generated by a horizontal transverse force {SH wave system). 
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check the agreement between the two computation techniques. Figures [5] and [6] show that the 
seismograms are essentially sensitive to porosity (f) in the uniform medium considered. The 
seismograms are also strongly influenced by perturbations of the consolidation parameter 
Cs, mineral density ps and shear modulus Gg- On the other hand, changes in fluid density 
Pj, mineral modulus Kg, fluid modulus Kf and permeability have only a weak influence 
on the wave amplitudes. We also note that the influence of parameters Cs and Gs on the one 
hand, and ps and p/ on the other hand, are very similar for the various force-displacement 
pairs. 

For the model with the lowest value of the consolidation parameter (corresponding to 




CC=92.0 % 



CC=99.3 % 



CC=99.8 % 



0.1 0.3 0.5 0.7 

Time (s) 

FIG. 4. Detailed comparison of individual seismic traces obtained with the Frechet derivative 
and with the discrete perturbation methods for various perturbation amplitudes. The model 
used for the computations is the 16-layer model presented in flgure [21 The seismic source 
is a vertical point force. The flgure shows the perturbation of the vertical displacement 
with respect to the consolidation parameter cs for a source-receiver offset of 500 m. The 
seismogram labeled a corresponds to the Frechet derivative formulation. The traces labeled 
6, c and d are respectively associated with discrete perturbations of the model parameter 
cs with relative amplitude contrasts of 1 %, 10 % and 50 %. The similarity between the 
Frechet derivative and discrete perturbation seismograms is quantifled by their correlation 
coefficient which is indicated above the traces. 
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the most consolidated material, figure [5]), the solid modulus Kg has a stronger influence on 
the seismograms than the fluid modulus Kf. On the contrary, for an unconsolidated medium 
(flgure[6]), the seismograms are mainly influenced by the fluid properties. With the model 
parameters used in this study, we flnd that the transition between these two regimes occurs 
for a consolidation parameter of 35. We also verifled, as suggested by our observations 
in section IV. B, that the porosity parameter shows the same behavior: the seismogram 
amplitudes in a high porosity medium depend more strongly on the fluid modulus than 




Fr-U Fr-V Fz-U Fz-V Ft-T 

force-displacement 

FIG. 5. Normalized maximum amplitude of the Frechet derivative seismograms for all solid 
force and displacement pairs. F^, Fr and Ft respectively denote the vertical, horizontal radial 
and horizontal tangential components of the forces acting on an average volume of porous 
medium, f/, V and T stand for the vertical, horizontal radial and horizontal tangential 
components of the solid displacement. The model parameters are those given in table 1 
except for the consolidation parameter Cg which is equal to 20. The amplitude scale is non 
linear for sake of readability. 
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on the solid modulus, and vice-versa for a medium with low porosity. As a consequence 
(and confirmation of field observations), the P- waves are strongly influenced by the fluid 
properties when they propagate in a fluid-saturated and poorly consolidated medium with 
high porosity. This makes it possible to determine the fluid characteristics from the seismic 
waveforms if these favorable conditions are met. Conversely, the estimation of the fluid 
characteristics will be more difficult in consolidated or unsaturated or low porosity media. 
Another consequence of the results shown in figures [5] and [6] is that and Cg are the most 
attractive parameters to invert for in an inversion procedure, assuming that ps and Gs can 
be estimated independently. 

B. Amplitude of the perturbation seismograms versus angle of incidence 

The previous figures already stressed the interdependence (or coupling) of some model 
parameters. Parameter coupling means that small perturbations of two or more parameters 




Fr-U Fr-V Fz-U Fz-V Ft-T 

force-displacement 



FIG. 6. Same as figure [5] with consolidation parameter Cg equal to 100. 
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FIG. 7. Energy of plane waves reflected from perturbations in p^, p/, /cq, 0, Kg, Kf, Gs 
and cs, as a function of incidence angle. The eight upper panels and eight lower panels 
respectively correspond to PP and 5*5* reflections. The curves are normalized with respect 
to the maximum value indicated above each panel. 
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result in similar modifications of the seismic response. An obvious consequence of param- 
eter coupling is that it becomes diffic ult or even impos sible to reliably estimate the model 
parameters in an inversion procedure (ITarantolal . 119861 ) . 



To look into this problem, we computed the plane wave responses corresponding to 10% 
perturbations of the model properties in the infinite medium described in table 1. Figure [7] 
shows the reflected energy for the eight model parameters as a function of angle of incidence 
at the perturbed layer, both for the PP and SS reflections. The smooth aspect of the curves 
is due to the intrinsic attenuation of the seismic waves propagating in the porous medium. 
The peaks and troughs of the curves are explained by th e strong variation s of th e reflection 



de la Cruz et al. 



mm . Some of 



and transmission coefficients, as shown for instance by 
these rapid variations are seen on the seismograms of figure [1] which were computed with 
the same model and source-receiver configuration. 



The magnitudes of the seismic responses shown in figure [7] are consistent with the study 
presented in section V.A. SS reflections are about twice as large as PP reflections for 
perturbations of the solid and fluid densities. For SS reflections, the maximum value is 
reached for a perturbation in porosity 0, whereas G^, Cg and ps produce perturbations of the 
same magnitude. For PP reflections, the most influential parameters are the solid density 
Ps and porosity 0. 



We also see in figure [7] that the radiation patterns associated with perturbations in p<j, 
Pf and fco on the one hand, and Kg and Kf on the other hand are exactly the same for the 
PP reflections. The same resemblance is observed for the SS reflections for the Ps, Pf and 
ko group of parameters, and for the Gg, and Cg group of parameters. We note in particular 
that permeability ko and densities pg and pf behave similarly despite their different roles 
in the constitutive equations. In all cases, the backscattered energy is maximum at normal 
incidence because of the shorter wave path and corresponding minimal attenuation. 
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VI. CONCLUSIONS 



We derived the Prechet derivatives of the seismic response of a depth-dependent porous 
medium. The Prechet derivatives are analytically expressed in terms of the Green's func- 
tions of the propagation medium through a perturbation analysis of the poro-elastic wave 
equations expressed in the plane wave domain. Started with a primary set of seven model 
parameters chosen because of their linear relationship with the wave equations, the deriva- 
tion was carried on with a secondary set of eight model parameters more convenient to use as 
physical parameters of the problem. The eight model parameters considered in our analysis 
are related to the fluid properties (density, bulk modulus), to the mineral properties (den- 
sity, bulk modulus, shear modulus) and to the arrangement of the porous material (porosity, 
permeabihty and consohdation parameter) . In the P — SV case, we derived Prechet deriva- 
tives for 3 different sources (horizontal and vertical point forces and explosive point source) , 
4 displacement components and 8 model parameters. In the SH case, we obtained 12 ex- 
pressions for 1 horizontal point force, 2 displacement components and 6 model parameters 
leading to non-zero Prechet derivatives. 

We checked the accuracy of these sensitivity operators in the time-distance domain by com- 
paring the waveforms computed from the first-order expressions with seismograms obtained 
by introducing discrete perturbations in the medium properties. The numerical tests were 
carried out both in a homogeneous and in a more complex earth model excited by oriented 
point forces. By and large, we found that our analytical expressions of the Prechet derivatives 
are remarkably accurate as long as the Born approximation assumptions are satisfied, that 
is, as long as the perturbations of the model parameters are weak and localized. However, as 
in other studies relying on the Born approximation, we showed that the first-order operators 
are robust enough to model parameter perturbations up to 20 % and layer thicknesses up 
to one fifth of the dominant wavelength. Purthermore, our formulation appears to be stable 
at all source-receiver offsets. 

Due to their analytical formulation, the sensitivity operators derived in this paper will be es- 
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pecially useful in full waveform inversion algorithms implemented with gradients techniques. 
As a first step toward such an application, we evaluated the sensitivity of the seismic re- 
sponse of a poro-elastic medium with respect to each model parameter. We showed that 
the porosity and consolidation parameter are the most attractive parameters to invert for, 
whereas the permeability appears to be the most difficult parameter to determine. The 
wave fields are more sensitive to the fluid bulk modulus than to the mineral bulk modulus, 
or inversely, according to porosity and consolidation parameter values. A multi-parameter 
inversion of backscattered energy looks challenging because of the strong coupling of several 
model parameters in a wide range of angles of incidence. Finally, this sensitivity study 
should prove useful for the interpretation of time-lapse monitoring surveys and for checking 
solutions (yet to come) accounting for 3D heterogeneities in the propagation medium. 
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